Die Dateiformate zur Speicherung von LiDAR-Daten: .las/.laz

LAS (Lidar Application Schema)

LAZ (LASzip Compressed Lidar)

Das R Paket lidR

“lidR ist ein R-Paket zur Manipulation und Visualisierung von Airborne Laser Scanning (ALS) Daten mit dem Schwerpunkt auf forstwirtschaftliche Anwendungen” (Roussel et al. 2022).

Import einer LAZ-Datei

las_files <- list.files(input_dir,
                        pattern = glob2rx("*.laz"),
                        full.names = TRUE)
  
las <- readLAS(las_files)

print(las)
## class        : LAS (v1.2 format 1)
## memory       : 2 Gb 
## extent       : 569000, 570000, 5708000, 5709000 (xmin, xmax, ymin, ymax)
## coord. ref.  : NA 
## area         : 1 kunits²
## points       : 36.47 million points
## density      : 36.47 points/units²
## density      : 15.01 pulses/units²

Validierung von LiDAR-Daten

las_check(las)

Plotting

plot(las)

max(las@data[["Z"]])
## [1] 741.042

Filterung der Punktwolke

  • nur Punkte kleiner gleich 500 m –> also könnten immer noch ca. 70 m hohe Bäume enthalten sein
las_new <- lidR::filter_poi(las, Z <= 500)
plot(las_new)

Cross section plot

p1 <- c(569250, 5708000)
p2 <- c(569750, 5709000 )
las_tr <- lidR::clip_transect(las_new, p1, p2, width = 4, xz = TRUE)

ggplot(las_tr@data, aes(X,Z, color = Z)) + 
  geom_point(size = 0.5) + 
  coord_equal() + 
  theme_minimal() +
  scale_color_gradientn(colours = height.colors(50))

Generierung eines Vegetationshöhenmodells (VHM)

Höhen Normalisierung

  • Beseitigung der Geländehöhe von der Gesamthöhe
  • zwei Ansätze:

1. Subtrahieren der DTM-Höhe von der DOM-Höhe

\(VHM = DOM - DGM\)

dtm <- lidR::rasterize_terrain(las_new, res = 1, knnidw())
plot(dtm, col = gray(1:50/50))

nlas <- las_new - dtm

hist(filter_ground(nlas)$Z, breaks = seq(-2.5, 2.5, 0.01), main = "", xlab = "Höhe")

  • nicht alle Bodenpunkte haben den Wert 0
  • Grund: Rasterformat des DGM –> Lage der Pixel stimmt nicht mit der Lage der Bodenpunkte überein
  • zu jedem Pixel mit einem Wert gibt es eigentlich mehrere Bodenpunkte
  • alle Punkte innerhalb eines bestimmten Pixels werden mit dem exakt gleichen Höhenwert normalisiert

2. Interpolation der Bodenpunkte und anschließend Subtraktion von den Nicht-Bodenpunkten

  • Interpolation der Höhe jedes einzelnen Punktes anhand von Bodenpunkten
nlas <- lidR::normalize_height(las_new, knnidw())

hist(filter_ground(nlas)$Z, breaks = seq(-0.5, 0.5, 0.01), main = "", xlab = "Höhe")

boxplot(nlas[["Z"]])

q1 <- quantile(nlas[["Z"]], .25)
q3 <- quantile(nlas[["Z"]], .75)
iqr <- IQR(nlas[["Z"]])

# Entfernung von Werten die das 1,5-fache des Interquartilsabstands größer als das dritte Quartil
# oder das 1,5-fache des Interquartilsabstands kleiner als das erste Quartil (Q1) betragen
nlas_rm_outliers <- subset(nlas, nlas[["Z"]] > (q1 - 1.5 * iqr) & nlas[["Z"]] < (q3 + 1.5 * iqr))
nlas_new <- subset(nlas_rm_outliers, nlas_rm_outliers[["Z"]] >= 0)

boxplot(nlas_new[["Z"]])

Rastern der normalisierten Punktwolke

  • verschiedene Methoden

Point-to-raster

  • Erstellung eines Gitters mit einer benutzerdefinierten Auflösung
  • Zuordnung der Höhe des höchsten Punktes zu jedem Pixel
chm_p2r <- lidR::rasterize_canopy(nlas_new, res = 1, algorithm = p2r())
plot(chm_p2r, col = height.colors(25))

  • Vorteil: rechnerisch einfach und schnell
  • Nachteil: einige Pixel können leer sein, wenn die Auflösung zu fein für die verfügbare Punktdichte ist
chm_p2r <- lidR::rasterize_canopy(nlas_new, res = 0.5, algorithm = p2r())
plot(chm_p2r, col = height.colors(25))

Lösung:

  • Verringerung der Anzahl leerer Pixel durch Ersetzung der Punkte mit Scheiben (mit definierten Radius)
  • Interpolation der verbleibenden leeren Pixel
chm_p2r <- lidR::rasterize_canopy(nlas_new, res = 0.5, algorithm = p2r(0.2, na.fill = tin()))
plot(chm_p2r, col = height.colors(25))

Triangulation

  • Erstellung eines unregelmäßiges Dreiecksnetzes (engl. Triangulated Irregular Network, TIN) mit den first returns
  • Interpolation innerhalb jeden Dreiecks, um einen Höhenwert für jedes Pixel eines Rasters zu berechnen
chm_tin <- lidR::rasterize_canopy(nlas_new, res = 0.5, algorithm = dsmtin())
plot(chm_tin, col = height.colors(25))

  • Vorteil: parameterfrei, keine leeren Pixel unabhängig von der Auflösung
  • Nachteil: Lücken und anderen Störungen (pits) durch first returns, die tief unter dem Kronendach liegen, geringere Performance bei Fehlen vieler Punkte

Lösung:

  • Nachbearbeitung des Höhenmodells (Glättung)
  • Entfernung von Dreiecken mit einer definierten Kantenlänge
chm_tin <- lidR::rasterize_canopy(nlas_new, res = 0.5, algorithm = dsmtin(max_edge = 8))
plot(chm_tin, col = height.colors(25))

Pit-free algorithm

  • aufeinanderfolgende Höhenschwellenwerte mit Delaunay-Triangulationen der first returns
  • Entfernung von zu großen Dreiecken für jeden Schwellenwert
  • Stapelung der einzelnen Raster und Beibehaltung der jeweils höchsten Pixel

chm_pitfree <- lidR::rasterize_canopy(nlas_new, res = 0.5,
                                      algorithm = pitfree(thresholds = c(0, 10, 20), max_edge = c(0, 1.5)))
plot(chm_pitfree, col = height.colors(25))

  • Vorteil: Vermeidung von pits, keine Nachbearbeitung oder Korrektur nötig
  • Nachteil: rechnerisch aufwändiger und langsamer

Vergleich der unterschiedlich berechneten VHMs

Quellen

Roussel J, Goodbody TR, Tompalski P (2022). The lidR package. https://r-lidar.github.io/lidRbook/index.html

Roussel J, Auty D, Coops NC, Tompalski P, Goodbody TR, Meador AS, Bourdon J, de Boissieu F, Achim A (2020). “lidR: An R package for analysis of Airborne Laser Scanning (ALS) data.” Remote Sensing of Environment, 251, 112061. ISSN 0034-4257, doi:10.1016/j.rse.2020.112061, https://www.sciencedirect.com/science/article/pii/S0034425720304314.